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Abstract. In this paper, we replace the asymptotic enrichments around the crack tip in the ex- 
^ | tended finite element method (XFEM) with the semi-analytical solution obtained by the scaled 
Q ■ boundary finite element method (SBFEM). The proposed method does not require special nu- 
merical integration technique to compute the stiffness matrix and it improves the capability of 
■ the XFEM to model cracks in homogeneous and/or heterogeneous materials without a priori 
'. knowledge of the asymptotic solutions. A heaviside enrichment is used to represent the jump 
<^ \ across the discontinuity surface. We call the method as the extended scaled boundary finite 
element method (xSBFEM). Numerical results presented for a few benchmark problems in the 
c* t ■ context of linear elastic fracture mechanics show that the proposed method yields accurate 
results with improved condition number. A simple MATLAB code is annexed to compute the 
terms in the stiffness matrix, which can easily be integrated in any existing FEM/XFEM code. 
Keywords: scaled boundary finite element method, extended finite element method, boundary integration, 
singular fields, partition of unity methods 
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Understanding the phenomenon of material failure is the key to design new materials to meet the increasing 
demand for high strength-to-weight and high stiffness-to-weight materials. In the classical approach, the 
governing equations are represented by partial differential equations (PDEs). Closed form or analytical 
CN solutions to these equations are not obtainable for most problems. Numerical methods such as the finite 
element method (FEM) [1], the boundary element method (BEM), meshless methods, spectral methods [2], 
discrete element methods [3], scaled boundary finite element method (SBFEM) [4, 5] are essential to 
analyze the equations. Amongst these, the FEM is a very popular and widely used approach, because of 
its versatility. The FEM relies on a conforming topological map and special elements have been developed 
for problems involving strong discontinuities. A constant remeshing is inevitable when the discontinuities 
evolve in time. It has been noted that the modification (augmentation, enrichment) of the finite element 
spaces will improve the behaviour. An intensive research over the past 3-4 decades has led to some of the 
robust methods available to model crack or discontinuities, for example extended (XFEM)/ Generalized 
FEM (GFEM)/ Partition of Unity FEM (PUFEM) [6, 7], hp- cloud [8] to name a few. The additional 
functions, also called the 'enrichment functions' are augmented through the partition of unity framework. 
These additional functions carry with them the information of the local nature of the problem. Since its 
inception, the XFEM has been applied successfully to moving boundary problems [9, 10, 11], to name a 
few. Although the XFEM is robust, it still suffers from certain difficulties and continuous research is carried 
out to improve the capability of the method. See Section that briefly discusses the various difficulties and 
attempts by various researchers to overcome the difficulties. For more detailed discussion, interested readers 
are referred to the literature [12] and references therein. 

Independently, Wolf and Song [5] developed a fundamental solution-less method, called the SBFEM 
for elasto-statics and elasto-dynamic problems. The SBFEM has emerged as an attractive alternate to 
model problems with singularities [13]. It combines the best features of the boundary element method and 
the finite element method. Numerical solutions are sought on the boundary, whilst the solution along the 



radial lines emanating from the 'scaling center' is represented analytically. Moreover, by utilizing the special 
features of the scaling center, the method allows the computation of stress intensity factor directly. 

Approach In this paper, we propose the xSBFEM, which combines the best features of the XFEM and 
the SBFEM, with the aim of improving the capability of the XFEM in representing the singular fields. As 
opposed to enriching a small region in the vicinity of the crack tip with asymptotic functions (known a 
priori either analytically or computed numerically), in the proposed method, that region is replaced with the 
semi-analytical solution obtained by the scaled boundary finite element method (SBFEM). The efficiency 
and the accuracy of the proposed method is illustrated by solving benchmark problems taken from linear 
elastic fracture mechanics. The effect of the size of the SBFEM domain on the global convergence in the 
energy is also studied numerically. Some of the salient features of the proposed method are: 

• Does not require a priori knowledge of the asymptotic expansions of the displacements and/or the 
stress fields. 

• As it does not require asymptotic enrichments, the need for special numerical integration technique 
is suppressed. 

• As only boundary information is required, the total number of unknowns is reduced. 

• The scaled condition number of the stiffness matrix of the xSBFEM is comparable and/or better 
than the XFEM. 

Outline The paper is organized as follows. In the next section, we briefly recall the basic equations 
of the XFEM. Section discusses the scaled boundary finite element method (SBFEM). The proposed 
technique, the eXtended SBFEM (xSBFEM) is presented in Section . Two different methods for computing 
the stress intensity factors, other than the interaction integral are discussed in Section . The efficiency and 
convergence properties of the proposed method are illustrated in Section with a few benchmark problems 
taken from linear elastic fracture mechanics, followed by some concluding remarks in the last section. 

Overview of the extended finite element method 

Consider SI C M 2 , the reference configuration of a cracked linearly isotropic elastic body (see Figure (1)). 
The boundary of f2 is denoted by T, is partitioned into three parts T U: T n and T c , where Dirichlet condition 
is prescribed on F u and Neumann condition is prescribed on T n and T c . The governing equilibrium equations 
for a 2D elasticity problem with internal boundary, r c defined in the domain Q and bounded by T is 

vJ<r + b = o g n (i) 

where V s (-) is the symmetric part of the gradient operator, is a null vector, cr is the stress tensor and b 
is the body force. The boundary conditions for this problem are: 
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where u = (u x ,u y ) T is the prescribed displacement vector on the essential boundary F u ; t = (t x ,t y ) T 
is the prescribed traction vector on the natural boundary T t and n is the outward normal vector. In this 
study, it is assumed that the displacements remain small and the strain-displacement relation is given by 
e = V s u. Let us assume a linear elastic behaviour, the constitutive relation is given by cr = D: e, where 
D is a fourth order elasticity tensor. Equation (1) is called the 'strong form' and computational solutions 
of Equation (1) rely on a process called 'discretization' that converts the problem into system of algebraic 
equations. The first step in transforming Equation (1) into a discrete problem is to reformulate Equation 



(1) into a suitable variational equation. This is done by multiplying Equation (1) with a test function. Let 
us define, 



U = {u£H 1 (tt) such that u| r „=u} 

V = {v£ff 1 (!l) such that v| r „ = 0} (3) 

The space U in which we seek the solution is referred to as the 'trial' space and the space V is called the 
'test' space. Let us define a bilinear form o(-, •) and a linear form £(■), 

o(u,v): = J cr(u): e(v) dfl; £(v) = J v • t dr (4) 
n r t 

where H 1 ^) is a Hilbert space with weak derivatives up to order 1. The weak formulation of the static 
problem is then given by: 

find u h eU h such that Vv h e V h a(u h ,v h ) = £(v h ) (5) 

where U h C U and V h C V. 




Fig. 1: Two-dimensional elastic body with a crack. 



extended Finite Element Method 

In the extended finite element method (XFEM), the classical FEM polynomial approximation space is 
extended by a set of functions, called the 'enrichment functions' . This is achieved through the framework of 
the partition of unity. By incorporating additional functions, the enriched FE basis can model singularities 
or internal boundaries [7]. 



Displacement approximation The displacement approximation can be decomposed into the standard 
part u^ td and into an enriched part u^ fcm as: 

u ft (x)=u\ td (x) + u^ m (x) (6) 

In the conventional XFEM, the nonsmooth behaviour around the crack tip is represented by enhancing 
the conventional finite element space by appropriate functions. For the case of linear elastic fracture me- 
chanics, two sets of functions are used: a Heaviside jump function to capture the jump across the crack 
faces and a set of asymptotic branch functions that span the two-dimensional asymptotic crack tip fields. 
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Fig. 2: A typical FE mesh with an arbitrary crack. 'Circled' nodes are enriched with d and 
'squared' nodes are enriched with near-tip asymptotic fields.' Reproducing elements' are the 
elements whose all the nodes are enriched. 



A generic form of the displacement approximation to capture the displacement jump across the crack faces 
and to represent the singular field is given by: 



u ft (x) = Y, ^(x)q/+ Y #j(*Mx)aj+ £ N K (x) ]T H a (x)b 



a=l 



u h std (x) 



(7) 



where jV c is the set of nodes whose shape function support is cut by the crack interior ('circled' nodes in 
Figure (2)) and AA f is the set of nodes whose shape function support is cut by the crack tip ('squared' 
nodes in Figure (2)). $ and E a are the enrichment functions chosen to capture the displacement jump 
across the crack surface and the singularity at the crack tip, q/ are the standard degrees of freedom, aj 
and bjj> are the nodal degrees of freedom corresponding to functions $ and H a , respectively. For example, 
a heaviside function is used to capture the jump in the displacement, whilst H Q is chosen to represent the 
near tip asymptotic fields. In the case of linear elastic fracture mechanics, these functions are given by: 
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For orthotropic materials, the asymptotic functions given by Equation (8) have to be modified because 
the material property is a function of material orientation. Asadpoure and Mohammadi [14, 15] proposed 
special near-tip functions for orthotropic materials as: 

{H a }i< a < 4 (r,6») = Vr <^ cos y ^/ 5 i(6»), cos y y/g 2 {0),sm-±-^/gi(O),sm-ly/g 2 (9) \ (9) 



where (r, 0) are the crack tip polar coordinates. The functions gi{i = 1,2) and 9i(i = 1,2) are given by: 



9j (6) = \^ C o S 2 9 + — —J , j = 1,2 

0j = arctan ^E^j . j = 1,2 (10) 

where ej(j = 1,2) are related to material constants, which depend on the orientation of the material [14]. 
The local enrichment strategy introduces the following four types of elements apart from the standard 
elements: 



• Split elements are elements completely cut by the crack. Their nodes are enriched with the discon- 
tinuous function i?. 

• Tip elements either contain the tip or are within a fixed distance independent of the mesh size. All 
nodes belonging to a tip element are enriched with the near-tip asymptotic fields. 

• Tip-blending elements are elements neighbouring tip elements. They are such that some of their 
nodes are enriched with the near-tip and others are not enriched at all. 

• Split-blending elements are elements neighbouring split elements. They are such that some of their 
nodes are enriched with the strongly or weakly discontinuous function and others are not enriched 
at all. 



Discretized Form By following Bubnov-Galerkin procedure, for arbitrary v h , we can write the dis- 
cretized form as: 
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(11) 



where K. uu ,~K aa and K ua ,K au are the stiffness matrix associated with the standard FE approximation, 
the enriched approximation and the coupling between the standard FE approximation and the enriched 
approximation, respectively. 



Remark The modification of the displacement approximation given by Equation (7) does not intro- 
duce a new form of the discretized finite element equilibrium equations, but leads to an enlarged problem 
to solve. 

Difficulties in the XFEM Although XFEM is robust and applied to a wide variety of moving boundary 
problems and interface problems, the flexibility provided by this class of methods also leads to associated 
difficulties: 



• Singular and discontinuous integrands When the approximation is discontinuous or non- 
polynomial in an element, special care must be taken to numerically integrate over enriched elements. 
One potential solution for numerical integration is to partition the elements into subcells (triangles 
for example) aligned to the discontinuous surface in which the integrands are continuous and dif- 
ferentiate [7, 16]. The purpose of sub-dividing into triangles is solely for the purpose of numerical 
integration and does not introduce new degrees of freedom. The other possible alternatives include: 
(a) Polar integration [17, 18]; (b) Complex mapping [19, 20, 21]; (c) Smoothed XFEM [22, 23, 24]; 
(d) Generalized quadrature [25, 26, 27]; (e) Equivalent polynomials [28] and (f) Adaptive integration 
techniques [29, 30, 31]. 



• Blending the different partitions of unity The local enrichment used in the conventional 
element leads to oscillations in the results over the elements that are partially enriched. This patho- 
logical behaviour has been studied in great detail in the literature. Some of the proposed techniques 
are: Assumed strain blending elements [32, 33], Direct coupling of the enriched and the standard 
regions [18, 32, 34], Corrected or weighted XFEM and Fast integration [35, 36], Hybrid crack ele- 
ment [37, 38], Hierarchical or spectral functions [39, 40]. 

• Ill-conditioning The addition of enrichment functions to the FE approximation basis could result 
in a severely ill-conditioned stiffness matrix. The consequence of this is that the stiffness matrix can 
become ill-conditioned [41, 42, 43, 44]. Some of the approaches proposed in the literature include 
Cholesky decompsition [41, 42, 17] of the submatrices of the stiffness matrix, domain decomposition 
based preconditioner [44], perturbation of the stiffness matrix [29, 45] and stable GFEM/XFEM [43, 
46]. 

• Additional unknowns With extrinsic enrichment, additional degrees of freedom (dofs) are in- 
troduced and the number of additional dofs depends on the number of enrichment functions and 
the number of such enrichments required. In case of the extrinsic enrichment, the approximation 
introduces additional unknowns (for example aj, b^-, see Equation (7)). These additional unknowns 
may increase the computational effort. By employing MLS technique on overlapping subdomains 
with ramp function [47] or by using crack tip elements developed for phantom node method [48, 49], 
the total number of additional unkowns can be reduced. 

It can be seen from the literature that different approaches have been proposed to overcome the afore- 
mentioned difficulties. In addition, Rethore eta/., [50] proposed to use analytical form of the displacement 
field in the vicinity of the crack tip, whilst retaining the conventional form in the rest of the domain. Rethore 
et al., used two overlapping domains to accomplish the task and the coupling between the two overlapping 
domains is established by matching the energy. The other approaches involve the work of Chahine and 
co-workers, viz., Cut-off XFEM [51] and spider XFEM [52]. Xiao and Karihaloo [38] employed the hybrid 
crack element with XFEM to improve the accuracy of the XFEM. The crack tip region was replaced with 
hybrid crack element and the XFEM was used to model the crack faces with jump functions. All of the 
above approaches, including the conventional XFEM, requires a priori knowledge of the asymptotic fields 
in an analytical form. Recently, Hattori et al., [53] derived new set of enrichment functions for anisotropic 
materials. To circumvent this problem, Menk and Bordas [44] proposed a numerical procedure to compute 
these enrichment functions, especially for anisotropic polycrystals. 

Remark Aforementioned difficulties are an artifact due to the asymptotic enrichments. In addition, 
the requirement for a priori knowledge of the asymptotic fields renders the application of the XFEM 
directly to heterogeneous materials for which the asymptotic fields do not exist in closed form or are 
very complex. Moreover, in some approaches, the computation of the stress intensity factors (SIFs) require 
further post-processing. Also, the introduction of the asymptotic fields, increases the complexity by requiring 
sophisticated numerical integration scheme to integrate over the enriched elements. 

Basics of the scaled boundary finite element method 

The scaled boundary finite element method (SBFEM) is a novel computational method developed by Wolf 
and Song [5], which reduces the governing partial differential equations to a set of ordinary differential 
equations. The SBFEM is semi-analytical and is suitable for solving linear elliptic, parabolic and hyperbolic 
partial differential equations. The SBFEM relies on transforming the conventional coordinate system to the 
scaled boundary coordinate system. Numerical solutions are sought around the circumferential direction 
using conventional finite element method, whilst in the radial direction smooth analytical solutions are 



obtained. Like the FEM, no fundamental solution is required and like the BEM, the spatial dimension is 
reduced by one, since only the boundary need to be discretized, resulting in a decrease in the total degrees 
of freedom. 

In the SBFEM, a scaling centre O is selected at a point from which the whole boundary of the domain 
is visible (scaling requirement). This requirement can always be satisfied by sub-structuring, i.e. dividing 
the structure into smaller subdomains. When modeling a cracked structure, a subdomain surrounding the 
crack tip is selected and the scaling centre is placed at the crack tip. The boundary of the subdomain is 
divided in to line elements (see Figure (3)). 




Fig. 3: A cracked domain modelled by SBFEM and the definition of local coordinate system, 
where the 'black' dots represent the nodes. 

The nodal coordinates on boundary are denoted as x&. As in a standard ID isoparametric finite element, 
the geometry of the element described by the coordinates x^r/), is expressed as 

MV) = N(r?)x fe (12) 

where N(?y) are the shape functions. Without loss of generality, the origin of the Cartesian coordinate 
system is chosen at the scaling centre. The geometry of the subdomain, described by x, is formed by 
scaling the boundary (Equation 12) 

x = fr b (v) (13) 

where £ is the normalized radial coordinate running from the scaling center towards the boundary, with 
£ = at the scaling center and £ = 1 on the boundary. The coordinates £ and rj are the so-called scaled 
boundary coordinates. They are related to the polar coordinates r and 0. The transformation is expressed 
as 



0(rj) = arctan 44 (14) 

where r& is the distance from the scaling centre to a point on the boundary. The transformation between 
the Cartesian coordinates and the scaled boundary coordinates is similar to the coordinate transformation 
in constructing isoparametric finite elements. 

Displacement approximation The displacements at any point is approximated by: 



u(e,f7) = N(f7)u(£) 



(15) 



where N(ry) are the shape functions of elements on the boundary and u(£) is the displacement along the 
radial lines, represented by a set of N analytical functions. By substituting Equation (15) in the definition 
of strain-displacement relations, the strains e(£,rj) are expressed as: 



where L is a linear operator matrix formulated in the scaled boundary coordinates 

L = bi (77)^ + ^2(77) 



(16) 
(17) 



and 
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(18) 



(19) 



(20) 



det(J) 

The determinant of the Jacobian matrix is: 

det(J) = x b {r})y h {r})^ - yb(v)xb(v),v 
where Xb(rj) and 7/6(77) are given by Equation (12). The stresses <r(^,r/) are then given by: 

<r(e,»7) = DBi(f7)u(e) j{ + |DB 2 (t ? )u(0 

where in the above equation, the definition of strain and the linear operator matrix given by Equation (16) 
- (17) are used with 

Bifa) = bi(j7)N(j7) 

B 2 fa)=b 2 fa)Nfa),„ (21) 

Upon substituting the Equations (20) and (16) in the virtual work statement and following the derivation 
in [5, 4], the virtual work yields 

5uJ ((E o eu(04 + E?u(0)| € =i-F) 

i . , 

f (22) 
- J 5u(0 T (E o e 2 u(0,^ + (E„ + Bf - Ei)£u(£),* - E 2 u(0) d£ = 

o 

where F is the equivalent boundary nodal forces and u& is the nodal displacement vector. By considering 
the arbitrariness of 5u(£), the following ODE is obtained: 

E o £ 2 u(0,« + (E + Bj - E!)eu(0, 5 - E 2 u(0 = (23) 
where E G ,Ei and E 2 are known as the coefficient matrices and are given by: 



E 
Ei 
E 2 



Bi( ?? ) T DBi(r / )det(J) d V , 
B 2 (7 1 ) T BB 1 (r ] )det(J) dr/, 
B 2 { V ) T T>B 2 (r])det(J) dr?. 



The boundary nodal forces are related to the displacement functions by 

F sbfcm = (E o eu(e), c + ET u (0)k=i 



(24) 
(25) 



Remark For a 2 noded line element, the coefficient matrices given by Equation (24) can be written 
explicitly. This is given in Appendix A. A MATLAB routine (EleCoeffMatrices2NodeEle.m) is given in 
Appendix B to compute the coefficient matrices. 

Computation of the stiffness matrix Equation (23) is a homogeneous second-order ordinary differ- 
ential equation, the general solution is given by [5, 4]: 



u(£) = cf>C x c = £cit 



(26) 



where cf) i and Aj can be interpreted as the deformation modes and the corresponding scaling factors that 
closely satisfy internal equilibrium in the £ direction [5, 4], and the constants Cj are the contribution of each 
mode. The deformation modes are obtained from the following quadratic eigenvalue problem obtained by 
substituting Equation (26) into Equation (23) as: 

(AfE o -Ai(E^-Ei)-E 2 )^ = (27) 

This equation yields 2N modes with pairs of eigenvalues (A, -A). Only the N modes with negative real 
parts of eigenvalue, which leads to finite displacements at the scaling centre, are selected. On the boundary, 
the displacement = u(£ = 1) is given by (Equation (26)): 

u b = 4>c (28) 

which permits the integration constants c to be determined from displacement on the boundary 

c = -1 u 6 (29) 

The nodal force on the boundary is obtained by substituting Equation (26) into Equation (25) as: 

F abfcm = (-E <£A + Ej<t>) c (30) 

Substituting Equation (29) into Equation (30) yields: 

_ sbfcm „sbfem . 

K u b = F (31) 
where |£ sbf * m j s t ne stiffness matrix given by: 

K Bbfom = -Eo^Xcf)' 1 + E^ (32) 



Alternatively, the quadratic eigen-problem in Equation (27) is rewritten as a standard eigen-problem of 
the matrix 

E^BJ — E~ 1 

El 1/ TT 1 — lTpT TTi XT' — 1 

2 + -tq-EJ Cj 1 — HilHi 

The eigenvalues and eigenvectors are partitioned as 
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(34) 



where A consists of all the eigenvalues with negative real parts. The displacement solution is expressed as 
(Equation (26)) 

u(£) = v n r A c (35) 
The boundary nodal force (Equation (30)) is written as 



The stiffness matrix is obtained as 



K 



v 2 ic 



v 2 iv n 



(36) 
(37) 



Remark A MATLAB routine (getSBFEMStiffMat.m) is given in Appendix B to compute the stiffness 
matrix, the deformation modes and the corresponding scaling factors. It is based on Equations (33) - (37). 

extended Scaled Boundary Finite Element Method 

The main idea of the proposed xSBFEM is to approximate the nonsmooth behaviour around the crack tip 
by a semi-analytical approach described in the previous section. In the proposed approach, the domain of 
Q is partitioned into three non-overlapping regions (Q = Q ™ U fi* ° m U iT ° m ) (see Figure (4)): 

fern 

• Q - contains a set of elements that are not intersected by the discontinuity surface. In this region 
the conventional FE approximation is employed to represent the smooth behaviour. 

• fT icm - contains a set of elements that are completely cut by the crack. Also called as the split 
element in XFEM terminology. The support of their nodes are enriched with a Heaviside function. 

• ff bi ° m - set of elements that define the SBFEM domain. The nonsmooth behaviour around the crack 
tip is modelled by a semi-analytical approach. 




0-/6A/-' 

Ale jV Bbfcm 



d enriched element (split element) 

Partially enriched element (blending element) 

Standard element 

Fig. 4: A typical FE mesh with a crack and different regions identified to represent the behaviour 
of the problem. In this descriptive figure, one layer of elements around the crack tip is replaced 
with the SBFEM domain. 'Circled' nodes are enriched with d and 'triangled' nodes are the 
nodes on the boundary of the Q* . Note that at the crack mouth, SBFEM requires two 
nodes, one on either side of the crack, whereas the discontinuity is captured by the enrichment 
function in case of the XFEM. 



Remark The SBFEM domain is constructed from the underlying FE mesh. After identifying a region 
around the crack tip, the information on the boundary of this region alone is required. The degrees of 
freedom of the unused nodes in the FE mesh are taken care of during the solution process. 

Salient features of the proposed approach The proposed approach aims at improving the accuracy of 
the XFEM in modeling problems with strong discontinuity. It shares the idea of representing the nonsmooth 
behaviour in the vicinity of the crack tip by analytical form similar to the work of Rethore et ai, [50] and 
Chahine et al., [51] but differs in basic approach and implementation. Some of the salient features are listed 
below: 

• An underlying FE mesh is used to define the SBFEM domain. 

• No a priori knowledge of the asymptotic fields is required. 

• In the vicinity of the crack tip, the integration is performed along the boundary of the SBFEM 
domain. 

• No special technique is required to match the different partitions of unity. 

• The existing approach can easily be incorporated into any existing FEM/XFEM code. 

Displacement approximation and the stiffness matrix The displacement approximation can be 
decomposed into the standard part u^ td and into an enriched part u^ fcm as given by Equation (7) without 
asymptotic enrichments. In the proposed approach, the standard part is further decomposed into the 
standard FEM Uf cm and into a semi-analytical part u^ bfcm . The displacement approximation is then given 
by: 



u h (x) 



£ N I (x)q J + £ JVj(x)i?(x)aj xe!]'~ur 
£ N K ( V )u bK (0 x = /(e,7 ? )Efi Bbf ° m 

^ eA f sbfcm 



where jV" fcm is the set of all nodes in the FE mesh that belongs to f2 em , jV sbfem j s the set of nodes in the 
FE mesh that belongs to Q" and jV" xfcm is the set of nodes that are enriched with the Heaviside function 
H. Nj and Nj are the standard finite element functions, q/ and aj are the standard and the enriched 
nodal variables associated with node / and node J, respectively, u^x are the nodal variables associated 
with the nodes on the boundary of $T ° m and Nk are the standard FE shape functions defined on the 
scaled boundary coordinates in Q° (see Section ). 

For $7 = f2 ™ U fT the discretized is given by Equation (11). For f2 = fT , the discretized form 
is given by Equation (31) and the stiffness matrix is given by Equation (32). Note that u^, in Equation 
(31) are the nodes on the boundary of SBFEM domain fT ° m that share with the boundary of and 
fT™. These nodes are identified from the underlying FE mesh. Although, the examples shown here are 
for structured mesh, the proposed technique is not limited to a structured mesh. Next, we describe the 
coupling of the three regions. 

Blending different partitions of unity Before discussing the coupling of different regions, let us define 
the notations used for the displacement variables in different regions. 

fern fcm sbfcm 

• q/ - standard degrees of freedom in O and u q s are the standard degrees of freedom in nfi 

• a/ - enriched degrees of freedom in $T' cm corresponding to the enrichment function $ that describes 
the jump in the displacement across the discontinuity surface. 



• MhK - degrees of freedom in fi sl "™ . This is further decomposed into two parts: are the degrees 
of freedom for the nodes in Q em n ff bem and u% K are the degrees of freedom for the nodes in 

,_sbfem _xfcm 

o nn . 

• u x p - contains the standard and the enriched degrees of freedom for the nodes in f] x '° m that share 
with the s ° m . 

Coupling of f2' cm , rT iom and ff bfem For the nodes that on the boundary of f2 f ° m and 0'"™, no special 
coupling technique is required. As the nodal unknown coefficients are required to be continuous across the 
boundary, the unknown coefficients from the SBFEM and the FEM are the same and are assembled in the 

fern xfcm 

usual way. A similar procedure is followed when assembling the U and ST to the global stiffness and 
to the global force vector. 

Coupling of Q Bbem and fT tom In case of the SBFEM, a crack in a body is modelled by an open domain 
as shown in Figure (3) with two nodes, one on either side of the crack, whilst, by using a heaviside function, 
a crack is modeled in the XFEM. The addition of a heaviside function introduces additional degrees of 
freedom that accounts for the jump in the displacement across the crack surface. Figure (5) shows a typical 
approach to model crack using the XFEM and the SBFEM. The displacement approximation for the XFEM 
and the SBFEM is given by: 



^■xfem 

i=i i=i 

UsVem = X>jfa)uw(0 (39) 



1=1 



where nf, is the total number of nodes on the boundary of ff As the displacements are required to be 
continuous, the displacements from the XFEM domain are matched to the displacements from the SBFEM 
domain by using a transformation matrix (see below). The continuity of the displacements are enforced 
only on the side that is common to both the ST"" and ff 



Numerical illustration The coupling between these two regions is illustrated below for a simple case 
where the crack is straight and the ff is to the right of fT"™ (see Figure (5)). The displacement 
approximation for the XFEM and the SBFEM is given by Equation (39). As explained earlier, to ensure 
compatibility of displacements, the displacements from the SBFEM domain are transformed to the dis- 
placements in the XFEM domain. This is done on the side that is common to both the regions, viz., sides 
connecting 2-3 in the XFEM domain and B-A: F-E in the SBFEM domain (see Figure (5)). From the XFEM 
approximation for the displacements, the displacements at the mouth of the crack, ie., at points l F' and 
'A' can be written as: 



u F (x, y) = N r(xF)<U + Yl N rM [ r ( x f) ~ a/ 

1=2,3 1=2,3 
= iV/(x F )q/ + 2N 2 (x F )a 2 

1=2,3 

u A (x, y) = J2 iV KxA)q/ + N i(*a) [H(x a ) ~ #(xj)] a 7 

1=2,3 1=2,3 

= N i(*a)cli - 2N 3 ( XA )a 3 (40) 
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Fig. 5: Coupling the XFEM and the SBFEM domain. 1 circled' nodes are the enriched nodes. 
1 Dashed' line represent the crack. 

Note that in Equation (40), the displacements up and are expressed in terms of nodal displacements 
q/ and the enriched displacements &j. Also note that these displacement should match with the displace- 
ments from the SBFEM domain. A transformation matrix can be derived from the above equation that 
translates the SBFEM nodal unknown displacements (viz, u^, ug, ue, up) to the XFEM nodal unknown 
displacements (viz, q 2 , q3, a 2 , a 3 ). In matrix form, 







I 











U E 







I 








U/l 


h 


iV 2 (x A ) 


iV 3 (x A ) 





-2iV 3 (x A ) 


k U F t 




N 2 (x F ) 


iV 3 (x F ) 


2iV 2 (x F ) 










1 













u bK 



Tu^ F 



(41) 



where I is the identify matrix and T is the transformation matrix that relates the displacements on the side 
that is common to ff and $T "\ To ensure the compatibility of the displacements and to assemble the 
stiffness matrix to the global stiffness matrix, the stiffness matrix in 0° and the force vector is rearranged 
as: 

K aa K a b 
Kf, a K hb 
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r < K 1 


H 


r F a 1 




l <K J 




1 F * i 



(42) 



where K aa , and F a corresponds to the stiffness matrix, the displacement vector and the force vector 
for the nodes that does not share with the element that is completely cut by the discontinuity surface in 
the J7 xcm , K^, uf K and F^ are the stiffness matrix, the displacement vector and the force vector that 
share with the element that is completely cut by the discontinuity surface and K a (,,K& a are the coupling 
matrices. Now, applying the transformation matrix, we get 
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(43) 



Remark The transformation matrix T requires only the evaluation of shape functions at the crack 
mouth, corresponding to the edge that is common to the fi"™ and ff ° m . 



Calculation of the stress intensity factors 



SIF from semi-analytical solution of displacements The relation between the stress intensity factors 
(SIFs) and the crack tip displacement fields for an isotropic material at a coordinate system as shown in 
Figure (3). The displacement field corresponding to the singular stress terms (superscript s for singular) is 
given by: 
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(44) 



where Kj and K// are the mode-l and mode-ll SIFs, G is the shear modulus, u is the Poisson's ratio and 
k is 

3 — 4f for plate strain 

(3 — i/)/(l + ^) for plane stress 

Consider the crack mouth opening displacement with r = r Q and 9 = ±7r, we have 
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(46) 



The singular stress terms can be identified by examining the eigenvalues of the terms in the semi- 
analytical solution (Equation (26)) obtained in the scaled boundary finite element method. For the case 
of a crack in a homogeneous subdomain, the eigenvalues of the singular terms are approximately equal to 
0.5. The displacement functions of the singular term is expressed as 



u s (£)=£ 



0.5 ^ , 

i=L II 



(47) 



The corresponding nodal displacement on the boundary (£ = 1) is written as 



u,, 
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(48) 



The values of the displacement of the singular stress term at the two nodes on the crack mouth are 
extracted from u| and substituted into the last term in Equation (46) leading to the stress intensity factors 
Kj and Kij. 

SIF from semi-analytical solution of stresses The stress intensity factors Kj and Kjj can be 

extracted from their definitions based on singular stresses a s yy and a s xy at the crack front 9 = 



Kj 
Kjj 



'2irr 



a yy (9 = 0) 

<y(0 = 0) 



(49) 



The singular stress modes are evaluated by using Equation (20) and the corresponding displacement 
modes in Equation (47) 



D[A i B 1 (? ? )+B 2 (? ? )]^ 



(50) 



As in the stress recovery in the finite element method, the computation is performed element-by-element 
at the Gauss integration point of an order lower by one than that for full integration. This yields an array of 



values of the singular stress modes at discrete values of angular coordinate 9 (Equation (14)). The values 
at the crack front 9 = are obtained by interpolation and denoted as *f>i(8 = 0). Along the crack front, 
£ = r/L a (Equation (14)) applies and the singular stresses are expressed as 



a(r,9 = 0) 



£ = 0) 

i=i,n 



(51) 



Substituting the stress components a yy and a xy obtained from Equation (51) into Equation (49) yields 
the stress intensity factors 



Remark The stress singularity is handled analytically and the standard stress recovery techniques 
are employed in extracting the stress intensity factors. 

Numerical Examples 

In this section, we illustrate the effectiveness of the proposed method by solving a few benchmark problems 
taken from linear elastic fracture mechanics. We first consider an infinite plate under tension, then a 
plate with an edge crack loaded with far field tension and shear. Then as a last example, we consider 
an orthotropic plate with a central crack and an edge crack with far field tension. In this study, unless 
otherwise mentioned, bilinear quadrilateral elements are used. The results from the proposed method are 
compared with those of the standard XFEM and the corrected XFEM. 

Infinite plate under remote tension In the first example, a plate with a crack loaded in remote 
tension is considered. Along the boundary ABCD (see Figure (6)), the closed form near tip displacements 
are imposed, given by (for plane strain conditions, see Equation (44)): 



where Kj = o^/tui denotes the stress intensity factor, a is the remote tension, v is the Poisson's ratio and 
E is the Young's modulus. All simulations are performed with a = 10 4 N/m 2 on a square with sides of 
length 10m. 

The scaled condition number & of the transformed stiffness matrix (K = P T KP) with increasing total 
degrees of freedom is shown in Figure (7). The matrix P is a preconditioner and it is chosen such that 
the condition number of the transformed matrix K is smaller than the original matrix. In this study, the 
diagonal scaling, also called as the Jacobi pre-conditioner is chosen, given by 



It can be seen from Figure (7) that the scaled condition number ^ increases with increasing degrees 
of freedom in case of all the approaches, viz., corrected XFEM (where blending correction with fixed area 
of enrichment is employed), standard XFEM ( 'topological' enrichment without blending correction) and 
proposed approach with different number of layers of SBFEM. But the rate of increase of the scaled 
condition number in case of the corrected XFEM is greater than the standard XFEM and the proposed 
approach. The scaled condition number of the proposed approach is lower than the standard XFEM and 
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(53) 



P = ^diag(K) 



(54) 




Fig. 6: Infinite plate with a center crack under remote tension. Region A — B — C — D is chosen 
as the computational domain and analytical displacements are prescribed along A — B — C — D 
given by Equation (53). 
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Fig. 7: Scaled condition number evaluated for different mesh densities. Std. XFEM refers to the 
conventional XFEM with 'topological' enrichment and the Corrected XFEM refers to XFEM 
with 'geometric' enrichment with blending correction. 




Fig. 8: Infinite plate under far field tension: Convergence results in the SIF. The rate of con- 
vergence is also given in the figure. 

decreases with increasing the number of layers of the SBFEM domain. The rate of convergence of the 
relative error in the SIF is shown in Figure (8) along with conventional XFEM with only tip element 
enriched ('topological enrichment') and XFEM with two layers of elements around the tip enriched with 
asymptotic fields. In this example, the numerical SIF is computed using the interaction integral. It is seen 
that the relative error in the numerical SIF decreases with increasing mesh density and that the proposed 
approach yields more accurate results for same number of enrichment layers. Moreover, the relative error 
in the numerical SIF decreases with increasing number of layers of SBFEM domain around the crack tip 
for same number of elements (see Figure (9)). 
Edge crack in tension and shear loading 

Edge crack in tension Consider a plate with an edge crack loaded by tension a = 1 over the top 
and bottom edges. The geometry, loading and boundary conditions are show in Figure (10). The reference 
mode I SIF is given by 

Ki = F (-|) a^E (55) 
where a is the crack length, H is the plate width and F (jj) is an empirical function given by (for ^ < 0.6) 

F (f ) - 1.12 - 0,3! (f ) + 10,5 - 2,12 (£)' + 30,0 (f )' ,53, 

The convergence of mode I SIF with mesh size and with the number of layers of SBFEM domain for 
a plate with an edge crack is shown in Figure (11). It can be seen that decreasing the mesh parameter h, 
the relative error in the numerical SIF decreases and that the error further decreases with increasing the 
size of the SBFEM domain. A similar behaviour is observed in the previous example. It can be concluded 
that 4 or 5 layers around the crack tip can yield accurate results. Figure (12) shows the convergence of 
the SIF with mesh refinement. The numerical SIF is computed by various methods, viz., stress based and 
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Fig. 9: Infinite plate with a center crack - convergence of numerical SIF with mesh density for 
different number of layers around the crack tip. 
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Fig. 10: Plate with an edge crack under tension 
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Fig. 11: Plate with an edge crack in tension: convergence of SIF with increasing mesh density 
for various number of layers around the crack tip. 

displacement based method (see Section ) and by the interaction integral method. It can be seen that with 
mesh refinement, the numerical SIF approaches the empirical value for all the different approaches. In this 
case, 5 layers of elements around the crack tip are replaced by the SBFEM domain. 

Edge crack in shear In this example, consider a plate with an edge crack subjected to a shear load 
r = 1 N/m 2 as shown in Figure (13). The material properties are assumed to be: Young's modulus, E = 
3xl0 7 N/m 2 and Poisson's ratio v = 0.25. The reference SIFs are given as: Kj = 34 N-m~ 3 / 2 and Ku = 
4.55 l\l-m~ 3//2 . From Table 1, it is evident that with mesh refinement, the computed SIFs converge to the 
reference mixed mode SIFs. It is noted that in this example, 5 layers of elements around the crack tip are 
replaced with the SBFEM domain. 

Table 1: Normalized mode I and mode II stress intensity factor for a plate with an edge crack 
with shear loading, where Kj mp = 34 N-m _3//2 and K e ^ v = 4.55 N-m~ 3 / 2 . 



Mesh mode I SIF (Kj/K e j mp ) mode II SIF (K n /K e ™ p ) 



20x40 


0.99662 


0.99419 


0.99500 


0.99253 


0.99301 


0.99488 


30x60 


0.99871 


0.99635 


0.99756 


0.99349 


0.99398 


0.99587 


40x80 


0.99968 


0.99737 


0.99870 


0.99385 


0.99435 


0.99629 


50x100 


1.00022 


0.99794 


0.99933 


0.99402 


0.99453 


0.99651 


60x120 


1.00056 


0.99829 


0.99972 


0.99413 


0.99464 


0.99664 



1.005 



LL 

CO 

<u 
■u 
o 
£ 

<D 
N 

15 
£ 



0.995 



0.99 



Empirical 

— Stress based 
-0— Displacement based 
-0— Interaction Integral 




15 20 25 30 35 40 45 50 

number of elements along x-direction 



55 



60 



65 



Fig. 12: Plate with an edge crack loaded in tension convergence results in the numerical SIF. 
Numerical SIF is computed using different techniques, viz., stress based, interaction integral 
and displacement based (see Section ). 




Fig. 13: Plate with an edge crack under shear: geometry and boundary conditions 



Crack in an orthotropic body 



Center crack in an orthotropic body Consider a square plate (h/w =1) with a center-crack of 
length 2a = 0.4 under uniform far field tension along the two opposite sides (see Figure (14). In this example, 
the Poisson's ratio and the shear modulus are taken as: v\ 2 = 0.03 and G\ 2 = 6 GPa, respectively. The 
Young's moduli E\ and E 2 are computed from the following expressions [53]: 

E x = G 12 (<S> + 21/12 + I) 

E, = § (57) 

where $ is the material parameter defined by the ratio between the Young's moduli. Figure (15) shows the 
influence of the material parameter on the normalized mode I SIF (Ki/(a^Jim)). The results are compared 
with those obtained in Ref. [53]. It is emphasized that in Ref. [53], new set of enrichment functions are 
derived for anisotropic materials and as such no enrichment functions are required in the proposed technique. 
The SIF in case of the proposed method is computed using the stress definitions as outlined in Section 
. It can be seen that the results obtained with the xSBFEM are in very good agreement with the results 
available in the literature. The SIF at the crack tips A and B are also shown in Figure (15) and they are 
similar as expected. 
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Fig. 14: A square plate with a center crack under uniform far field tension. E\ and E 2 are the 
Young's moduli along the x— and the y— directions. 



Edge crack under uniform tension Consider a square plate (h/w = 1) with a double edge-crack 
(a/w = 0.5). The plate is subjected to uniform tractions, as shown in the Figure (16). The plate is 
a symmetric angle ply composite laminate consisting of four graphite-epoxy laminae, with the following 
elastic properties: E\ = 144.8 GPa, E 2 = 11.7 GPa, G 12 = 9.66 GPa and u 12 = 0.21. In this case, the 
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Fig. 15: Influence of material parameters on the normalized mode I SIF (Ki j '{a^pKa)) for an 
orthotropic plate with a center crack. 

influence of the fiber orientation on the SIF is studied by varying the fiber angle from 9 = 0° to 9 = 
90°. Due to symmetry, only one half of the plate is modelled. A structured quadrilateral mesh of 60x120 
elements is used. Figure (17) shows the variation of normalized mode I SIF {K i / {a ^/^w J )) with the fiber 
orientation. It can be seen that with increasing fiber orientation, the numerical SIF increases and reaches 
a maximum value at 9 =90°. The mode I SIF is symmetric with respect to the fiber orientation 9 =90° as 
expected. 



Conclusions 

In this paper, we have proposed a new method by combining the scaled boundary finite element method 
with the extended finite element method for problems with strong discontinuity. Instead of enriching a 
small region in the vicinity of the crack tip with asymptotic fields, the stiffness of the small region around 
the tip is computed by employing the SBFEM. With a few examples from linear elastic fracture mechanics, 
the effectiveness of the proposed method is illustrated. It is seen that the proposed method outperforms 
the conventional approach to handle discontinuities inside the domain. The key features of the proposed 
method without asymptotic enrichment are: (a) eliminates the need to sub-divide the element containing 
the crack tip; (b) suppresses the need to integrate singular functions usually appearing in the stiffness 
matrix, as there is no enrichment with asymptotic fields; (c) the entire process can be parallelized; (d) the 
scaled condition number of the proposed technique is better than the XFEM counterpart and (e) reduces 
the computational effort as the total number of unknowns are reduced. Although the numerical examples 
and the formulation presented here are for a structured mesh with bilinear elements. The proposed approach 
can easily be extended to other element formulations, viz., triangular, higher order quadrilateral or polygonal 
elements and can be applied to heterogeneous materials. 
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Fig. 17: Variation of the normalized mode I SIF (Ki/^a^/iw)) with orthotropic angle 
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Appendix A 

Computation of the coefficient matrices and the stiffness matrix For a 2-node line element with nodal 
coordinates (x\, y\), (x 2 , 2/2) an d the shape functions along the line element are given by: 

N 1 ( V ) = ^(l- V y } N 2 (r,) = ^(l + V ). 

The integrations over rj in the coefficient matrices E OJ Ei and E2 (see Equation (24)) are performed 
analytically, yielding 
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where the following abbreviations are introduced 
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Appendix B 

function [ Kb, lambda, v] = getSBFEMStiffMat ( coord, conn, D) 

% purpose: compute the stiffness matrix for the SBFEM domain 

% 

% Inputs: 

% coord — nodal coordinates 



% conn — element connectivity along the boundary 

% D — elasticity matrix 

% 

% Outputs 

% Kb — stiffness matrix 

% lambda, v— eigenvalues and eigenvectors 

% lambda, v are required later for computing SIFs 

% 

nd = numel ( coord ) ; id = l:nd; 

EO = zeros(nd, nd); El = zeros(nd, nd ) ; E2 = zeros(nd, nd ) ; 
for ie = 1 : s i ze ( conn , 2 ) 

x = coord ( 1 , con n ( 1 : 2 , i e ) ) ; y = coord (2 , conn ( 1 : 2 , i e ) ) ; 

d = [ 2* con n ( 1 , i e ) — 1 2*conn(l,ie) 2* con n ( 2 , i e ) — 1 2*conn(2,ie) ]; 

[ eO_ele , el.ele, e2_ele] = E I eCoeff M a t ri ces2 N od e E I e (x , y, D); 

EO(d,d)=EO(d,d)+eO_ele ; El(d , d)=El(d , d)+el_ele ; E2(d ,d)=E2(d ,d)+e2_ele ; 

end 

m=E0\[El' -eye(nd)]; Z= [m; El*m( : , id )-E2-ld-12*E0 El*m( : , nd+id ) ] ; 

[v, d] = eig(Z); lambda = diag(d); 
%sort eignvalues in ascending order 

[~ , idx] = so rt ( rea I ( lambda ),' ascend '); 
% rearrange eigenvalues and eigenvectors 

lambda = lambda ( idx ( id )) ; v = v(:, idx(id)); 
% stiffness matrix 

Kb = real (v(nd+id , :)/v(id , :)); 

function [ eO , el, e2] = E I eCoeff M a t r i ces2 N od e E I e ( x, y, D ) 

% purpose: compute elemental coefficient matrices. 

% 

% Inputs: 

% x,y— current element nodal coordinates 

% D — elasticity matrix 

% 

% Outputs : 

% eO , el , e2 — SBFEM element coefficient matrices 

% 

% 

a = x(l)*y(2)-x(2)*y(l); 

CI = [ y(2)-y(l) 0; -(x(2)-x ( 1 ) ) ; -(x(2)-x(l)) y(2)-y(l) ]; 
C2 = 0.5*[ y(2) + y(l) 0; -(x(2) + x ( 1 ) ) ; -(x(2) + x(l)) y(2) + y(l)]; 
Q0 = 0.25/a*(Cl'*D*Cl); Ql = -0.25/a *(C2 ' *D*C1 ) ; Q2 = . 25 / a * ( C2 ' * D*C2 ) ; 
eO = 2/3*[ 2*Q0 Q0 ; Q0 2*Q0] ; 

el = l/3*[ -Q0 Q0; Q0 -Q0 ] + 2*[ -Ql -Ql; Ql Ql ]; 
e2 = l/3*[ Q0 -Q0; -Q0 Q0 ] + 4*[ Q2 -Q2; -Q2 Q2 ]; 

References 

[1] Zienkiewicz OC, Taylor RL, Zhu JZ. The finite element method: its basics and fundamentals. 
6 th edn., Elsevier Butterworth Heinemann, 2000. 



[2] Boyd JP. Chebyshev and Fourier Spectral Methods. Dover Publications Inc, 2000. 



Munjiza AA. The combined finite-discrete element method. Wiley Publishers, 2004. 

Deeks A, Wolf J. A virtual work derivation of the scaled boundary finite-element method for 
elastostatics. Computational Mechanics 2002; 28:489-594. 

Wolf J, Song C. The scaled boundary finite-element method - a fundamental solution-less 
boundary element method. Computer Methods in Applied Mechanics and Engineering 2001; 
190:5551-5568. 

Babuska I, Melenk J. The partition of unity finite element method. International Journal for 
Numerical Methods in Engineering 1997; 40:727-758. 

Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. Inter- 
national Journal for Numerical Methods in Engineering 1999; 45:601-620. 

Garcia 0, Fancello E, Barcellos C, Duarte C. Hp clouds in mindlin's thick plate model. Inter- 
national Journal for Numerical Methods in Engineering 2000; 47:1381-1400. 

Duddu R, Bordas SPA, Chopp D, Moran B. A combined extended finite element and level set 
method for biofilm growth. International Journal for Numerical Methods in Engineering Jan 
2008; 74(5):848-870, doi:10.1002/nme.2200. 

Aubertin P, Rethore J, de Borst R. Dynamic crack propagation using a combined molecular 
dynamics/extended finite element approach. International Journal for Multiscale Computational 
Engineering 2010; 2:221-235. 

Chessa J, Smolinski P, Belytschko T. The extended finite element method (XFEM) for so- 
lidification problems. International Journal for Numerical Methods in Engineering Jan 2002; 
53(8):1959-1977, doi:10.1002/nme.386. 

Belytschko T, Gracie R, Ventura G. A review of extended /generalized finite element methods for 
material model. Modelling and simulation in Materials Science and Engineering 2009; 17:1-24. 

Song C, Tin-Loi F, Gao W. A definition and evaluation procedure of generalized stress intensity 
factors at cracks and multi-material wedges. Engineering Fracture Mechanics 2010; 77:2316- 
2336. 

Asadpoure A, Mohammadi S. A new approach to simulate the crack with the extended finite 
element method in orthotropic media. International Journal for Numerical Methods in Engi- 
neering 2007; 69:2150-2172. 

Asadpoure A, Mohammadi S. Crack analysis in orthotropic media using the extended finite 
element method. Thin Walled Structures 2007; 44:1031-1038. 

Bordas S, Nguyen VP, Dunant C, Guidoum A, Nguyen-Dang H. An extended finite element 
library. International Journal for Numerical Methods in Engineering 2007; 71:703-732. 

Bechet E, Minnebo H, Moes N, Burgardt B. Improved implementation and robustness study of 
the X-FEM for stress analysis around cracks. International Journal for Numerical Methods in 
Engineering 2005; 64(8):1033-1056, doi:10.1002/nme.l386. 



[18] Laborde P, Pommier J, Renard Y, Salaun M. High-order extended finite element method for 
cracked domains. International Journal for Numerical Methods in Engineering September 2005; 
64(3):354-381, doi:10.1002/nme.l370. 

[19] Natarajan S, Bordas S, Mahapatra DR. Numerical integration over arbitrary polygonal domains 
based on schwarzchristoffel conformal mapping. International Journal for Numerical Methods 
in Engineering 2009; 80:103-134. 

[20] Natarajan S, Mahapatra DR, Bordas S. Integrating strong and weak discontinuities without 
integration subcells and example applications in an xfem/gfem framework. International Journal 
for Numerical Methods in Engineering 2010; 83:269-294. 

[21] Natarajan S, Baiz P, Bordas S, Kerfriden P, Rabczuk T. Natural frequencies of cracked func- 
tionally graded material plates by the extended finite element method. Composite Structures 
2011; 93:3082-3092. 

[22] Bordas S, Natarajan S, Kerfriden P, Augarde C, Mahapatra DR, Rabczuk T, Pont SD. On 
the performance of strain smoothing for quadratic and enriched finite element approximations 
(XFEM/GFEM/PUFEM). International Journal for Numerical Methods in Engineering 2011; 
86:637-666. 

[23] Bordas S, Rabczuk T, Nguyen-Xuan H, Nguyen VP, Natarajan S, Bog T, Quan DM, Nguyen 
VH. Strain smoothing in FEM and XFEM. Computers & Structures 2010; 88:1419-1443. 

[24] Baiz PM, Natarajan S, Bordas S, Kerfriden P, Rabczuk T. Linear buckling analysis of cracked 
plates by SFEM and XFEM. Journal of Mechanics of Materials and Structure 2011; 6:1213- 
1238. 

[25] Mousavi SE, Sukumar N. Generalized Gaussian quadrature rules for discontinuities and crack 
singularities in the extended finite element method. Computer Methods in Applied Mechanical 
and Engineering 2010; 199:3237-3249. 

[26] Mousavi SE, Sukumar N. Generalized duffy transformation for integrating vertex singularities. 
Computational Mechanics 2010; 45:127-140. 

[27] Mousavi SE, Xiao H, Sukumar N. Generalized Gaussian quadrature rules on arbitrary polygons. 
International Journal for Numerical Methods in Engineering 2010; 82:99-113. 

[28] Ventura G. On the elimination of quadrature subcells for discontinuous functions in the ex- 
tended finite-element method. International Journal for Numerical Methods in Engineering 
2006; 66(5):767-795. 

[29] Strouboulis T, Babuska I, Copps K. The design and analysis of the generalized finite element 
method. Computer Methods in Applied Mechanical and Engineering 2000; 181:43-69. 

[30] Liu XY, Xiao QZ, Karihaloo BL. XFEM for direct evaluation of mixed mode SIFs in homo- 
geneous and bi-materials. International Journal for Numerical Methods in Engineering 2004; 
59:1103-1118. 

[31] Xiao QZ, Karihaloo BL. Improving the accuracy of XFEM crack tip fields using higher or- 
der quadrature and statically admissible stress recovery. International Journal for Numerical 
Methods in Engineering 2006; 66(9): 1378-1410. 



[32] Gracie R, Wang H, Belytschko T. Blending in the extended finite element method by discon- 
tinuous Galerkin and assumed strain methods. International Journal for Numerical Methods in 
Engineering 2008; 74:1645-1669. 

[33] Chessa J, Wang H, Belytschko T. On the construction of blending elements for local partition 
of unity enriched finite elements. International Journal for Numerical Methods in Engineering 
2003; 57:1015-1038. 

[34] Chahine E, Laborde P, Renard Y. A non-conformal eXtended Finite Eelement approach: Integral 
matching Xfem. Applied Numerical Mathematics 2011; 61:322-343. 

[35] Fries TP. A corrected XFEM approximation without problems in blending elements. Interna- 
tional Journal for Numerical Methods in Engineering 2008; 75:503-532. 

[36] Ventura G, Gracie R, Belytschko T. Fast integration and weight function blending in the ex- 
tended finite element method. International Journal for Numerical Methods in Engineering 
2009; 77:1-29, doi:10.1002/nme.2387. 

[37] Xiao QZ, Karihaloo BL. Direct evaluation of accurate coefficients of the linear elastic crack tip 
asymptotic field. Fatigue & Fracture of Engineering Materials & Structures 2003; 26:719-729. 

[38] Xiao QZ, Karihaloo BL. Implementation of hybrid crack element on a general finite element 
mesh and in combination with XFEM. Computer Methods in Applied Mechanical and Engi- 
neering 2007; 196:1864-1873. 

[39] Tranacon JE, Vercher A, Giner E, Feuenmayor FJ. Enhanced blending elemetns for XFEM 
applied to linear elastic fracture mechanics. International Journal for Numerical Methods in 
Engineering 2009; 77:126-148. 

[40] Legay A, Wang HW, Belytschko T. Strong and weak arbitrary discontinuities in spectral finite 
elements. International Journal for Numerical Methods in Engineering 2005; 64:991-1008. 

[41] Strang G, Fix GJ. An Analysis of the Finite Element Method. Prentice-Hall, 1973. 

[42] Fix GJ, Gulati S, Wakoff Gl. On the use of singular functions with finite element approximations. 
Journal of Computational Physics 1973; 13:209-228. 

[43] Babuska I, Banerjee U. Stable generalized finite element method (SGFEM). Technical Report 
11-07, Institute for Computational Engineering and Sciences 2011. 

[44] Menk A, Bordas S. A robust preconditioning technique for the extended finite element method. 
International Journal for Numerical Methods in Engineering 2011; 85:1609-1632. 

[45] Strouboulis T, Copps K, Babuska I. The generalized finite element method. Computer Methods 
in Applied Mechanical and Engineering 2001; 190:4081-4193. 

[46] Natarajan S, Banerjee U, Bordas S. extended finite element method: Implementation, accuracy 
and convergence properties. International Journal for Numerical Methods in Engineering 2012; 
Article in review. 

[47] Fries TP, Belytschko T. The intrinsic partition of unity method. Computational Mechanics 
2007; 40:803-814. 



[48] Song JH, Areias PM, Belytschko T. A method for dynamic crack and shear band propaga- 
tion with phantom nodes. International Journal for Numerical Methods in Engineering 2006; 
67:868-893. 

[49] Rabczuk T, Zi G, Gerstenberger A, Wall WA. A new crack tip element for the phantom- 
node method with arbitrary cohesive cracks. International Journal for Numerical Methods in 
Engineering 2008; 75:577-599. 

[50] Rethore J, Roux S, Hild F. Hybrid analytical and extended finite element method (HAX-FEM) 
: A new enrichment algorithm for cracked solids. International Journal for Numerical Methods 
in Engineering 2010; 81:269-285. 

[51] Chahine E, Laborde P, Renard Y. Crack tip enrichment in the (XFEM) using a cutoff function. 
International Journal for Numerical Methods in Engineering 2008; 75:629-646. 

[52] Chahine E, Laborde P, Renard Y. Spider-XFEM: an extended finite element variant for partially 
unknown crack-tip displacement. European Journal of Computational Mechanics 2008; 15:625- 
636. 

[53] Hattori G, Rojas-Diaz R, Saez A, Sukumar N, Garcia-Sanchz F. New anisotropic crack-tip 
enrichment functions for the extended finite element method. Computational Mechanics 2012; 
50:591-601. 



